/**
* @Arquivo: lagrangian.c
* 
*	Geração e Solução do problema de coloração usando a formulção de representantes 
*
* @author: Marcio Costa Santos
* @data: 20/05/2010
*
*/

#include "lagrangian.h"
#include "pl.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define TAG 2

extern int n;//o numero de vertices
extern int m;//o numero de arestas
extern int k;//o k da coloração
extern long va[];//a ordem dos verties
int nt;//numero de threads

CPXENVptr env[MAX_SIZE];//o ambiente do problema
CPXLPptr prob[MAX_SIZE];//os problemas
double x[MAX_SIZE][MAX_SIZE];//o vetor de pesos
double bestLB;//o melhor limite inferior encontrado
double bestUB;//o melhor limite superior encontrado
double lb;//o limite inferior
double ub;//o limite superior
double pi;//passo do lagrangeano
int noI;//melhoria?
int handle;//auxiliar
double coef[MAX_SIZE];//os coeficientes da função objetivo
long map[MAX_SIZE*MAX_SIZE]; //mapeia as variaveis nos seus mu

int rank;//rank do processo
int size;//tamanho da rede
MPI_Request res;//requição
MPI_Status sts;//status
int scheduling[MAX_SIZE];//escalonamento das tarefas

/*Função de ordenação para o qsort*/
int compar(const void* a, const void* b){
	int v = *(int*)a;
	int u = *(int*)b;
	return mu[v] - mu[u];
}

int define(){
	int i,status;
	pi = 2.0;
	noI = 0;
	lb = ub =  0;
	bestLB = (-1)*n;
	bestUB = n;
	for(i = 0; i < n; i++){
		mu[i] = Lx[i] = 2.0;
		env[i] = CPXopenCPLEX(&status);
	}
	if(status!= 0){
		printf("problemas para abrir cplex");
		fflush(NULL);
		exit(0);
		MPI_Finalize();
	}
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	MPI_Comm_size(MPI_COMM_WORLD, &nt);

	scheduling[0] = 0;
	for(i = 1; i < nt; i++){
		scheduling[i] = i*(n/(nt-1));
	}
	scheduling[nt-1] = n;

	return initPL();
}

int makePL(long v){
	return addConstrainsForVertex(v, &(env[v]) ,&(prob[v]));
}

double calculateUB(){
	int V[n];
	int cor[n];
	int i,j;
	int color = 1;
	int adj[n];
	int cand = -1;
	for(i = 0; i < n; i++){
		V[i] = i;
		cor[i] = 0;
		adj[n] = 0;
	}
	qsort(V, n, sizeof(int), compar);
	cor[0] = 0;
	for(i = 1; i< n; i++){
		for(j = 0; j < i; j++){
			if(!hasEdge(V[i], V[j])){
				adj[cor[j]] = 1;
			}
		}
		for(j = color - 1; j >= 0; j--){
			if(adj[j] == 0)cand = j;
		}

		if(cand == -1){
			cand = color;
			color++;
		}
		cor[i] = cand;

		cand = -1;
		for(j = 0; j < n; j++){
			adj[j] = 0;
		}

	}
	if(bestUB > k*color){
		bestUB = k*color;
	}
	return k*color;
}

int noImproviment(){
	noI++;
	if(noI > 30){
		pi = pi/2.0;
		noI = 0;
	}
	return 0;
}


double calculateLB(){
	int i;
	double value = 0.0;
	for(i = 0; i < n; i++){
		value += Lx[va[i]];
		value += k*mu[va[i]];
	}
	if(value > bestLB){
		bestLB = value;
	}	
	else noImproviment();
	return value;
}

int updateMuLambda(){
	double g[MAX_SIZE];
	double theta = 0.0;
	double theta_div = 0.0;
	int i,j;
	for(i = 0; i < n; i++){
		g[va[i]] = (double)k;
		for(j = 0; j <= i; j++){
			g[va[i]] -= k*x[va[j]][va[i]];
		}
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*((1.2)*ub - lb)/theta_div;
	else theta = pi/n;
	//calcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] += theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
			if(g[va[i]] <= 0.0) g[va[i]] = 0.0;
		}
	}
	theta_div = 0.0;
	//recalcula subgradente
	for(i = 0; i < n; i++){
		theta_div += g[va[i]]*g[va[i]];
	}
	if(theta_div >=  0.001)theta = pi*((1.05)*ub - lb)/theta_div;
	else theta = pi/n;
	//calcula subgradientes
	for(i = 0; i < n; i++){
		mu[va[i]] += theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
		}
	}
	return 0;
}


double solvePL(int v){
	int i = 0;
	Lx[va[v]] = k - k*solve(env[v], prob[v], v) - k*mu[va[v]];
	getX(env[v], prob[v], x[va[v]], v);
	if(Lx[va[v]] > 0 ){
		Lx[va[v]] = 0;
		for(i = 0; i < n; i++){
			x[va[v]][i] = 0.0;
		}
	}
	return Lx[va[v]];
}

int addObjFunc(long v){
	int i, j;
	j=0;
	for(i = v+1; i < n; i++){
		if(hasEdge(va[i], va[v])){
			coef[j++] = mu[va[i]];
		}
	}
	return changeObj(env[v],prob[v], coef, v);
}


double lagrangian(double* sol){
	int i;
	for(i = 0; i < n; i++){
		makePL(i);
	}
	if(rank == 0){
		while(pi > 0.005){
			for(i = 0; i < nt;i++){
				sendMu(i);
			}
			receiveSol();
			lb = calculateLB();
			ub = calculateUB();
			updateMuLambda();
		}
		endProcs();
	*sol = bestUB;
	}else{
		int active = 1;
		while(active){
			double data[n+1];
			receiveUpdate(data, &active);
			sendSol();
		}
		exit(0);
		MPI_Finalize();
	}

	return bestLB;
}

void receiveUpdate(double* data, int* active){
	MPI_Recv(data,n+1,MPI_DOUBLE,0,MPI_ANY_TAG,MPI_COMM_WORLD,&sts);
	if((int)data[n] == 0){
		int i;
		for(i = 0; i < n; i++){
			mu[i] = data[i];
		}
		int init = scheduling[rank-1];
		int end = scheduling[rank];
		for(i = init; i < end; i++){
			addObjFunc(i);
			solvePL(i);
		}
	}
	else{
		*active = 0;
	}
}

void sendMu(int i){
	int j;
	double data[n+1];
	for(j = 0; j < n; j++){
		data[j] = mu[j];
	}
	data[n] = 0;
	MPI_Send(data,n+1,MPI_DOUBLE,i,TAG,MPI_COMM_WORLD);
}


void receiveSol(){
	int i,j,l;
	double data[n+2];
	for(i = 1; i < nt; i++){
		int init = scheduling[i-1];
		int end = scheduling[i];
		for(j = init; j < end; j++){
			MPI_Recv(data,n+2,MPI_DOUBLE, i, TAG ,MPI_COMM_WORLD, &sts);
			int v = (int)data[n+1];
			for(l = 0; l <n; l++){
				x[va[v]][l] = data[l];
			}
			Lx[va[v]] = data[n];
		}
	}

}

void endProcs(){
	int data[n+1];
	data[n] = 1.0;
	MPI_Bcast(data,n+1,MPI_DOUBLE,0,MPI_COMM_WORLD);
}

void sendSol(){
	int init = scheduling[rank-1];
	int end = scheduling[rank];
	int i,j;
	double data[n+2];
	for(i = init; i < end; i++){
		for(j = 0; j < n; j++){
			data[j] = x[va[i]][j];
		}
		data[n+1] = i;
		data[n] = Lx[va[i]];
		MPI_Send(data,n+2,MPI_DOUBLE,0,TAG,MPI_COMM_WORLD);
	}

}

